
class ProbitTemp{
	public:
		ProbitTemp(double *Y,mat X, mat s){
			_Y=Y;
			_X=X;
			_s=s;
			_p=X.n_cols;
			_n=X.n_rows;
			_xi_n=100;
			_is=inv(s);
			d= new boost::math::normal_distribution<>(0,1); 
		}
		~ProbitTemp()
		{
			delete d;
		}
		double Likelihood(mat theta1,double xi)
		{
			mat theta=theta1.t();
			double L=0;
			double sum=0;
			int m=theta.n_cols;
			//Xb
			int n=Get_n();
			mat X=Get_X();
			double *Y=Get_Y();
			for(int i=0;i<n;i++)
			{	
			//	cout << theta.col(0) << "\\";
				mat foo= X.row(i);
				sum=dot(foo,theta)/xi;
			//	cout << sum << " ";
				double phi=Phi(sum);
				double phi2=1-Phi(sum);
				if(Y[i]==1){
					L+=log(phi);
				}else if(Y[i]==0){
					L+=log(phi2);
				}else{
					cout << "Warning values of Y!= {0,1}";
				}
			}
			//cout << exp(L);
			return L;
		}
		double Prior(mat theta)
		{
			double res=0;
			mat foo=theta.t();
			int d=theta.n_elem;
			mat bar=foo*_is*foo.t();
			//cout << "//" << bar;
			res=-d*0.5*(log(2*PI))-0.5*log(det(_s))-0.5*as_scalar(bar);
			return 0;

		}	
		inline double Density_utn(mat theta, int n){
			return Likelihood(theta,_xi_n)+Prior(theta);
		}
		inline double Weight(mat theta,double *W,int n)
		{
			double temp=0;
			double temp2=0;
			double temp3=0;
			//cout << "2 " << bar << " 3";
			int M=theta.n_cols;
			double *foo=new double[M];
			for(int i=0;i<M;i++)
			{	
				temp=(this->Likelihood(theta.col(i),_xi_n))-(this->Likelihood(theta.col(i),_xi_n1));
				W[i]*=exp(temp);	
				//cout << "\n";
			}
			delete[] foo;
		}
		inline double* Eta_star(double xi, double *W, mat theta, int Bool)
		{
			double temp=0;
			double temp0=0;
			double temp2=0;
			double temp3=0;
			double xin=_xi_n1-xi;//-_phi_n1;
			//cout << "2" << bar << "3";
			int M=theta.n_cols;
			double *foo=new double[M];
			if(Bool==1)
			{
			mat Temp(M,1);
			for(int i=0;i<M;i++)
			{	
				Temp(i,0)=(this->Likelihood(theta.col(i),xin))-(this->Likelihood(theta.col(i),_xi_n1));
				foo[i]=W[i]*exp(Temp(i,0));	
				/*temp0=(this->Likelihood(theta.col(i)));
				temp2=(this->Prior(theta.col(i)));
				temp3=-log(_Start->d(theta.col(i)));
				cout << "temp2 " << temp2 <<"\n";
				cout << "temp0 " << temp0 <<"\n";
				cout << "temp3 " << temp3 <<"\n";*/
				//cout << "foo: "<< foo[i] << " ";
			}
			_Targ=Temp;
			}else{
				for(int i=0;i<M;i++)
				{	
					foo[i]=W[i]*exp(_Targ(i,0));	
				}
			}
			//cout << foo <<"\n";
			return foo;
		}
		inline double Phi(double x)
		{
			double t=cdf((*d),x);
			if(t==1)
			{
				t=0.99999999999;
			}
			if(t==0)
			{
				t=0.00000000001;
			}
			return t;
			//return boost::math::cdf<double, boost::math::policies::policy<> >(boost::math::normal_distribution(0,1),x);
		}
		inline void Set_xi(double xin){
			_xi_n1=_xi_n;
			_xi_n=xin;
			
		}
		inline int Set_p(int p){ _p=p;}
		inline int Get_p(void){return _p;}
		inline double Get_bn(void){return _xi_n;}
		inline int Get_n(void){return _n;}		
		inline	mat Get_X(void){return _X;}
		inline	double* Get_Y(void){return _Y;}
		inline mat Get_v(void){return _s;}
	private:
		double _xi_n;
		double _xi_n1;
		mat _X;
		mat _Targ;
		double *_Y;
		boost::math::normal_distribution<> *d;
		int _p;
		int _b;
		mat _s;
		mat _is;
		int _n;
		Distribution::Distribution *_Start;
};
